Esta guía explica cómo realizar algunas tareas comunes de análisis de datos de recuento de secuenciación microbiana utilizando tidytacos. Ilustraremos utilizando un conjunto de datos con muestras de microbioma humano del tracto vaginal, tomadas de este artículo por Lebeer et al.

Tidytacos: la filosofía del paquete

Un objeto tidytacos es simplemente una lista de tres tablas:

El paquete se llama tidytacos porque cada una de las tablas está ordenada (en inglés tidy): cada fila representa una observación y cada columna una variable (puedes encontrar más información sobre la ordenación de datos en este enlace).

Configuración

En caso de que aún no hayas instalado tidytacos, puedes instalarlo usando devtools:

install.packages("devtools")
devtools::install_github("LebeerLab/tidytacos")

Para esta guía, sólo necesitamos cargar tres paquetes: tidytacos (por supuesto), el conjunto de paquetes tidyverse, y VTutorials.

library(tidyverse)
library(tidytacos)
library(VTutorials)

Exploremos vdata

Nuestro conjunto de datos de ejemplo está disponible en el paquete VTutorials y no es necesario importarlo ni convertirlo. Se llama “vdata” y comenzamos inspeccionando la tabla de muestras:

glimpse(vdata$samples)
## Rows: 200
## Columns: 28
## $ sample                               <chr> "SAMEA13411869", "SAMEA13410951",…
## $ sample_id                            <chr> "s16", "s100", "s114", "s170", "s…
## $ Technical.run_20201009_isala_cross_1 <lgl> FALSE, TRUE, FALSE, FALSE, TRUE, …
## $ Technical.run_20201016_isala_cross_2 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201023_isala_cross_3 <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, …
## $ Technical.run_20201030_isala_cross_4 <lgl> FALSE, FALSE, FALSE, TRUE, FALSE,…
## $ Technical.run_20201106_isala_cross_5 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201120_isala_cross_6 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201204_isala_cross_8 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201211_isala_cross_9 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20210119_isala_cross_7 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.libsize                    <dbl> 37957, 22782, 32150, 7707, 24123,…
## $ General.Age                          <dbl> 24, 60, 24, 33, 29, 24, 23, 42, 3…
## $ Health.BMI                           <dbl> 30.11938, 19.00391, 25.03414, 23.…
## $ Health.Antibiotic.3months            <lgl> FALSE, TRUE, TRUE, TRUE, FALSE, F…
## $ Sexual.Intercourse.24hours           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ plate                                <chr> "PlateBA", NA, "PlateBB", "PlateB…
## $ volume                               <dbl> 2.2025248, NA, 1.5580467, 1.83601…
## $ dna_conc                             <dbl> 11.169, NA, 15.789, 20.648, 24.65…
## $ read_conc                            <dbl> 17233.4038, NA, 20634.8110, 4197.…
## $ read_conc_norm                       <dbl> 0.9693078, NA, 1.1606229, 0.57860…
## $ lib_size                             <dbl> 37957, NA, 32150, 7707, 24123, 29…
## $ `Sample month`                       <dbl> 7, NA, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
## $ `Transit time`                       <dbl> 1, NA, 1, 1, 1, 2, 1, 1, 6, 2, 1,…
## $ valencia_subcst                      <chr> "I-A", "IV-C0", "III-A", "IV-B", …
## $ valencia_cst                         <chr> "I", "IV-C", "III", "IV-B", "IV-B…
## $ max_genus1                           <chr> "L. crispatus", "Prevotella", "L.…
## $ max_genus2                           <chr> "Gardnerella", "Other", "Prevotel…

Luego echamos un vistazo rápido al número total de muestras, ASV y lecturas (reads en inglés) en el objeto tidytacos:

tacosum(vdata)
## n_samples    n_taxa   n_reads 
##       200       942   5022096

Hagamos un diagrama de barras apilado de un subconjunto de muestras

Podemos crear muy fácilmente un gráfico para explorar un subconjunto de nuestras muestras (por ejemplo, solo muestras de participantes que han tomado antibióticos en los últimos 3 meses y son menores de 40 años) de la siguiente manera:

vdata %>%
  filter_samples(Health.Antibiotic.3months == TRUE, General.Age <= 40) %>%
  tacoplot_stack()

La función filter_samples hace lo que dice: filtrar muestras. También eliminará los taxones de la tabla de taxones que tengan cero lecturas totales en las muestras restantes. La función tacoplot_stack devuelve una buena visualización de diagramas de barras apiladas de los taxones más abundantes en nuestras muestras.

¿Se te ocurren otras ideas de como filtrar y rápidamente visualizar vdata?

Subconjunto de datos

Nuestra siguiente pregunta para este conjunto de datos es hasta qué punto la actividad sexual dentro de las 24 horas antes de tomada las muestras afecta la composición microbiana de las participantes en edad reproductiva (asumamos <=40 años). Para hacernos una idea primero lo podemos visualizar:

vdata_less40 <- vdata %>% filter_samples(General.Age <= 40)
tacoplot_stack(vdata_less40)+
  geom_point(aes(y=-0.02,color=Sexual.Intercourse.24hours))

## Diversidad Alpha Para explorar la diversidad alfa, creemos una versión enrarecida (rarefied en inglés) del conjunto de datos:

vdata_rar <- vdata %>%
  add_total_count() %>%
  filter_samples(total_count >= 2000) %>%
  rarefy(2000) %>%
  add_alpha()

La función add_total_count agregará números totales de lectura de muestra a la tabla de muestra.

La función rarefy submuestreará aleatoriamente todas las muestras n veces. Solo funciona si el recuento de lecturas de cada muestra es igual o superior a n. 

Para determinar la riqueza de ASV, optamos por enrarecer primero, pero esto puede depender de sus datos.

La función add_alpha se puede utilizar para agregar varias métricas de diversidad alfa a la tabla de muestra.

Podemos visualizar la diversidad alpha de muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no:

vdata_rar %>%
  samples() %>%
  ggplot(aes(x = Sexual.Intercourse.24hours, y = observed, fill = Sexual.Intercourse.24hours)) +
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(height = NULL)

Análisis de coordenadas principales (PCoA)

Nos gustaría abordar las diferencias entre muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no. También estamos más interesados en los géneros que en los ASV.

Una PCoA podría ofrecer información:

vdata_genus <- vdata %>%
  aggregate_taxa(rank = "genus")

tacoplot_ord_ly(vdata_genus, Sexual.Intercourse.24hours, samplenames = sample, 
                dim = 3)

La función aggregate_taxa fusiona todas las filas de la tabla de taxones en un nivel taxonómico específico, en este caso el nivel de género. Como ocurre con todas las funciones de tidytacos, todas las demás tablas del objeto tidytacos se ajustan en consecuencia.

La función tacoplot_ord_ly determinará la abundancia relativa de taxones en las muestras y luego utilizará las diferencias de Bray-Curtis para ordenar muestras en un espacio bidimensional (o tridimensional) según su composición taxonómica. La adición argumental de plotly “_ly” hace que la figura sea interactiva, lo cual es realmente bueno para el trabajo exploratorio. Esto también funciona para otras funciones de visualización.

Relación entre la composición de la comunidad microbiana y variables.

La siguiente pregunta lógica es hasta qué punto la actividad sexual reciente (dentro de las 24 horas) determina la variabilidad de la composición de la comunidad microbiana. No olvidemos que la composición microbiana vaginal varia con la edad, por eso incluyamos edad en el modelo.

perform_adonis(vdata_genus, c("General.Age", "Sexual.Intercourse.24hours"))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = as.formula(paste("counts_matrix", formula_RHS, sep = " ~ ")), data = metadata, permutations = permutations)
##                             Df SumOfSqs      R2      F Pr(>F)    
## General.Age                  1    1.717 0.02676 5.5016  0.001 ***
## Sexual.Intercourse.24hours   1    0.956 0.01489 3.0607  0.018 *  
## Residual                   197   61.500 0.95835                  
## Total                      199   64.173 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La función perform_adonis realizará un ANOVA PERMutacional para determinar el efecto de las variables de muestra en las diferencias de Bray-Curtis de las comunidades. El resultado muestra que la edad es un contribuyente a la composición de la comunidad microbiana (R cuadrado = 0.02676).

Análisis de abundancia diferencial

A continuación, nos gustaría saber cuáles de los 20 géneros más abundantes son significativamente más abundantes en las participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no.

vdata_genus <- vdata_genus %>% add_codifab(Sexual.Intercourse.24hours, 
                                           max_taxa = 20)
vdata_genus$taxon_pairs <- filter(vdata_genus$taxon_pairs, wilcox_p < 0.05)
tacoplot_codifab(vdata_genus, FALSE_vs_TRUE)

La función add_codifab agregará una tabla llamada taxon_pairs al objeto tidytacos, con la abundancia diferencial del taxón entre las dos condiciones (con respecto al taxón de referencia), para cada par de un taxón y un taxón de referencia.

La función tacoplot_codifab devuelve un gráfico para visualizar la abundancia diferencial de taxones entre condiciones, en comparación con todos los demás taxones como referencia. Podemos observar que es más probable que Staphiloccus y L. iners group sea típico encontrar cuando una participante ha tenido actividad sexual dentro de las 24 horas antes de tomada las muestras.

Es de destacar que existen muchos métodos de análisis de abundancia diferencial y ninguno de ellos es perfecto. Interprete tus resultados con cuidado.

Descargo de responsabilidad: Esta guía es una adaptación de la versión en inglés del tutorial de tidytacos.